Chapter 10 Dominance-microbiome analysis
10.0.1 Posterior estimates
The next step is to get the posterior estimates of beta, which is the parameter that links dominance with the microbiome. We employ a support threshold of ‘0.05’ to assign significance, meaning that parameters with posterior overlaps of <10% are considered significantly different.
# Select desired support threshold
support=0.9
negsupport=1-support
# Basal tree
postestimates_tree <- genome_tree %>%
keep.tip(., tip=m$spNames)
# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
mutate(value = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
pivot_wider(names_from = variable, values_from = value) %>%
rename(intercept=2) %>%
select(genome,dominance,groupvariable) %>%
column_to_rownames(var="genome")
# Aggregate basal GIFT into elements
function_table <- genome_gifts %>%
to.elements(., GIFT_db) %>%
as.data.frame()
#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% postestimates_tree$tip.label) %>%
arrange(match(genome, postestimates_tree$tip.label)) %>%
mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
column_to_rownames(var = "genome") %>%
select(phylum)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% postestimates_tree$tip.label) %>%
arrange(match(genome, postestimates_tree$tip.label)) %>%
select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
select(colors) %>%
pull()
# Basal ggtree
postestimates_tree <- postestimates_tree %>%
force.ultrametric(.,method="extend") %>%
ggtree(., size = 0.3)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
scale_fill_manual(values=colors_alphabetic)+
labs(fill="Phylum")
#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()
# Add posterior significant heatmap
postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.2, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
labs(fill="Trend")
#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()
#Add functions heatmap
postestimates_tree <- gheatmap(postestimates_tree, function_table, offset=0.6, width=3.5, colnames=FALSE) +
vexpand(.08) +
coord_cartesian(clip = "off") +
scale_fill_gradient(low = "#f4f4f4", high = "steelblue", na.value="white") +
labs(fill="Function fullness")
#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()
# Add completeness barplots
postestimates_tree <- postestimates_tree +
geom_fruit(data=genome_metadata,
geom=geom_bar,
grid.params=list(axis="x", text.size=2, nbreak = 1),
axis.params=list(vline=TRUE),
mapping = aes(x=length, y=genome, fill=completeness),
offset = 3.85,
orientation="y",
stat="identity") +
scale_fill_gradient(low = "#cf8888", high = "#a2cc87") +
labs(fill="Genome completeness")
postestimates_tree +
vexpand(.25, 1) # expand top
### Positively associated genomes
getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="dominance", trend=="Positive") %>%
arrange(-value) %>%
left_join(genome_metadata,by=join_by(genome==genome)) %>%
select(genome,phylum,class,order,species,value) %>%
tt()| genome | phylum | class | order | species | value |
|---|---|---|---|---|---|
| A_DRC05_bin.59 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__Duncaniella sp910577345 | 1.000 |
| C_HT_bin.102 | p__Bacillota | c__Bacilli | o__Lactobacillales | s__Limosilactobacillus reuteri | 0.998 |
| A_HTC12_bin.48 | p__Bacillota | c__Bacilli | o__Lactobacillales | s__Lactobacillus johnsonii | 0.995 |
| C_HT_bin.9 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__QXXE01 sp910589115 | 0.994 |
| A_DTC03_bin.77 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__CAG-873 sp910586975 | 0.992 |
| B_C04M3_bin.80 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__ | 0.991 |
| B_C05M2_bin.12 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__CAG-873 sp910578095 | 0.991 |
| A_CDC10_bin.30 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__Duncaniella sp910589485 | 0.990 |
| B_C03M5_bin.33 | p__Bacillota_B | c__Dehalobacteriia | o__UBA4068 | s__ | 0.990 |
| A_CDC13_bin.65 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Eubacterium_J sp009774535 | 0.987 |
| A_T3C12_bin.34 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__CAG-485 sp009775355 | 0.985 |
| C_DR_bin.56 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__CAG-485 sp910578075 | 0.985 |
| B_C04M3_bin.108 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Acutalibacter sp009936035 | 0.984 |
| B_C06M3_bin.61 | p__Bacillota | c__Bacilli | o__Erysipelotrichales | s__C-19 sp910576855 | 0.981 |
| B_C12F1_bin.72 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__Duncaniella sp910578515 | 0.980 |
| C_HR_bin.153 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.980 |
| A_T1C10_bin.3 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__UBA7173 sp002491305 | 0.979 |
| B_C11F5_bin.68 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__UBA7173 sp001689485 | 0.979 |
| B_C12F1_bin.61 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__CAG-873 sp009775535 | 0.977 |
| A_DRC04_bin.49 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Enterenecus sp910587255 | 0.975 |
| C_DR_bin.4 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Avispirillum sp011957885 | 0.974 |
| A_HTC03_bin.32 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Pelethomonas sp910587645 | 0.973 |
| B_C04M5_bin.101 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__Paramuribaculum sp910579675 | 0.973 |
| B_C11F3_bin.29 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__CAG-95 sp009911035 | 0.973 |
| A_CDC12_bin.65 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__Coproplasma sp003979335 | 0.972 |
| B_C13F3_bin.77 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__CAG-873 sp910587915 | 0.972 |
| A_CDC13_bin.13 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__UBA3402 sp002358555 | 0.971 |
| B_C05M3_bin.38 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Eubacterium_J sp910574915 | 0.971 |
| B_C13F3_bin.53 | p__Bacillota | c__Bacilli | o__Erysipelotrichales | s__Erysipelatoclostridium cocleatum | 0.970 |
| A_DTC13_bin.55 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__RACS-045 sp910579785 | 0.969 |
| A_T2C12_bin.11 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__NSJ-51 sp910585415 | 0.969 |
| B_C06M5_bin.62 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__CAG-475 sp910577815 | 0.968 |
| B_C13F3_bin.5 | p__Bacillota_A | c__Clostridia | o__TANB77 | s__ | 0.967 |
| B_C13F3_bin.97 | p__Bacillota_A | c__Clostridia | o__TANB77 | s__ | 0.967 |
| C_CR_bin.18 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__VSOB01 sp910589075 | 0.966 |
| C_DT_bin.159 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__ | 0.965 |
| A_T2C05_bin.3 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__Muribaculum gordoncarteri | 0.964 |
| B_C13F2_bin.43 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Enterenecus sp910587255 | 0.963 |
| B_C04M3_bin.128 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Anaerotruncus sp000403395 | 0.962 |
| B_C13F3_bin.35 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__MGBC136627 sp910585935 | 0.961 |
| A_OPC13_bin.50 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__UBA1405 sp910589145 | 0.960 |
| B_C05M1_bin.75 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Clostridium_Q sp009911305 | 0.959 |
| B_C13F3_bin.70 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.959 |
| A_CDC05_bin.89 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__ | 0.958 |
| A_OPC10_bin.79 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__UBA3402 sp910574345 | 0.957 |
| C_OP_bin.44 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Ventrimonas sp910586205 | 0.957 |
| C_HT_bin.137 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.956 |
| C_CR_bin.157 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__ | 0.954 |
| C_T3_bin.64 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__NSJ-51 sp910585415 | 0.954 |
| C_DR_bin.31 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Anaerotignum sp910576545 | 0.953 |
| A_OPC10_bin.40 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.952 |
| B_C12F2_bin.47 | p__Bacillota | c__Bacilli | o__Lactobacillales | s__Enterococcus_D casseliflavus | 0.952 |
| C_HR_bin.13 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Evtepia sp910584805 | 0.952 |
| C_T2_bin.55 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__ | 0.952 |
| B_C04M2_bin.33 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Enterenecus sp910585265 | 0.951 |
| B_C04M3_bin.31 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Angelakisella sp013316495 | 0.951 |
| A_OPC03_bin.64 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Acetatifactor sp910589655 | 0.950 |
| A_T1C11_bin.9 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__CAG-485 sp002362485 | 0.949 |
| B_C04M3_bin.18 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Eubacterium_J sp910586445 | 0.949 |
| B_C03M3_bin.127 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Acutalibacter sp009936035 | 0.948 |
| B_C04M4_bin.45 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.948 |
| A_T3C03_bin.17 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__QANA01 sp910588425 | 0.947 |
| B_C04M2_bin.76 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__CAG-56 sp004793585 | 0.945 |
| A_HRC11_bin.35 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.944 |
| B_C04M4_bin.5 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Choladocola sp009774145 | 0.944 |
| C_DT_bin.37 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Acutalibacter sp009936035 | 0.944 |
| C_DR_bin.21 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__MD308 sp910578455 | 0.943 |
| A_CDC03_bin.78 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Lawsonibacter sp910577805 | 0.942 |
| B_C13F5_bin.28 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__Duncaniella sp910576785 | 0.942 |
| A_HRC10_bin.28 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__Gallimonas sp910578065 | 0.941 |
| A_T3C11_bin.23 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__RGIG4284 sp910576205 | 0.941 |
| A_HTC12_bin.21 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__UBA3402 sp910575585 | 0.940 |
| A_T2C05_bin.28 | p__Bacillota | c__Bacilli | o__Erysipelotrichales | s__MGBC163490 sp910588075 | 0.940 |
| B_C10F2_bin.5 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.940 |
| C_HR_bin.146 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Choladocola sp910575445 | 0.940 |
| B_C05M5_bin.57 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Pelethomonas sp910587785 | 0.939 |
| B_C10F5_bin.45 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Dysosmobacter sp000403435 | 0.939 |
| B_C12F3_bin.66 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__ | 0.939 |
| B_C11F5_bin.87 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__ | 0.938 |
| A_CRC05_bin.32 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Eubacterium_J sp910579075 | 0.937 |
| B_C04M1_bin.34 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Pelethomonas sp910579965 | 0.937 |
| B_C04M3_bin.90 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__14-2 sp910580015 | 0.937 |
| C_T2_bin.30 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__ | 0.934 |
| C_T3_bin.41 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Acetatifactor sp910584235 | 0.934 |
| C_CD_bin.124 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Pelethomonas sp910579965 | 0.933 |
| C_OP_bin.73 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__UBA3402 sp910585435 | 0.933 |
| C_HR_bin.50 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__UBA3402 sp910588325 | 0.932 |
| A_CDC03_bin.64 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__ | 0.930 |
| B_C03M3_bin.121 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Lachnoclostridium_B sp910574185 | 0.930 |
| B_C10F1_bin.59 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Choladocola sp910583895 | 0.929 |
| B_C10F5_bin.36 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__UBA3282 sp003611805 | 0.928 |
| B_C11F5_bin.21 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__COE1 sp000403335 | 0.928 |
| A_CDC10_bin.13 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__MGBC113161 sp910587565 | 0.927 |
| A_DRC13_bin.41 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Acutalibacter muris | 0.927 |
| C_CD_bin.62 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__1XD42-69 sp009911505 | 0.925 |
| C_OP_bin.15 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Roseburia sp910584215 | 0.925 |
| A_CDC05_bin.91 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.923 |
| C_T2_bin.92 | p__Bacillota_A | c__Clostridia | o__TANB77 | s__ | 0.923 |
| A_CDC05_bin.30 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.922 |
| A_CRC04_bin.71 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Sporofaciens sp910574885 | 0.922 |
| B_C13F3_bin.99 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Merdisoma sp910574255 | 0.922 |
| B_C04M1_bin.8 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Merdisoma sp910576325 | 0.921 |
| B_C03M2_bin.44 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__Duncaniella muris | 0.920 |
| B_C06M5_bin.16 | p__Actinomycetota | c__Coriobacteriia | o__Coriobacteriales | s__Adlercreutzia muris | 0.920 |
| C_T3_bin.99 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Pelethomonas sp910577915 | 0.920 |
| B_C05M1_bin.6 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Ventrimonas sp910584845 | 0.919 |
| A_DRC06_bin.37 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__CAG-873 sp009775265 | 0.918 |
| B_C12F1_bin.38 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Dysosmobacter sp910588005 | 0.918 |
| B_C10F1_bin.10 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Ventrimonas sp910577765 | 0.917 |
| B_C13F3_bin.93 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Anaerotruncus sp003612625 | 0.917 |
| A_HRC05_bin.71 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__MGBC100320 sp910577995 | 0.915 |
| A_T3C13_bin.3 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__ | 0.914 |
| B_C03M1_bin.60 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__MGBC100798 sp910584975 | 0.911 |
| C_OP_bin.100 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__ | 0.910 |
| C_DT_bin.155 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__ | 0.909 |
| A_T1C03_bin.24 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__UMGS1370 sp910577645 | 0.907 |
| B_C12F1_bin.23 | p__Bacillota_A | c__Clostridia | o__Oscillospirales | s__Pelethomonas sp910579965 | 0.907 |
| A_T2C13_bin.14 | p__Bacteroidota | c__Bacteroidia | o__Bacteroidales | s__UBA7173 sp001689685 | 0.906 |
| A_CRC05_bin.73 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__MGBC100320 sp910588305 | 0.905 |
| B_C03M1_bin.29 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Ventrimonas sp009911065 | 0.905 |
| B_C04M5_bin.14 | p__Actinomycetota | c__Coriobacteriia | o__Coriobacteriales | s__Adlercreutzia caecimuris | 0.905 |
| C_HR_bin.28 | p__Bacillota_A | c__Clostridia | o__Lachnospirales | s__Sporofaciens sp910585725 | 0.903 |
| B_C04M5_bin.93 | p__Bacillota_A | c__Clostridia | o__Christensenellales | s__Gallimonas sp910588465 | 0.902 |
10.0.2 Negatively associated genomes
getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="dominance", trend=="Negative") %>%
arrange(value) %>%
left_join(genome_metadata,by=join_by(genome==genome)) %>%
select(genome,phylum,class,order,species,value) %>%
tt()| genome | phylum | class | order | species | value |
|---|---|---|---|---|---|
| A_T3C05_bin.41 | p__Pseudomonadota | c__Gammaproteobacteria | o__Burkholderiales | s__Parasutterella excrementihominis | 0.095 |
10.1 Predict responses to dominance
Based on the fitted model, we can predict the structure of the microbial community for each level of dominance for a certain value of sequencing depth.
# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")
gradient = seq(0,1,0.1)
gradientlength = length(gradient)
#Treatment-specific gradient predictions
pred <- constructGradient(m,
focalVariable = "dominance",
non.focalVariables = list(logseqdepth=list(1)),
ngrid=gradientlength) %>%
predict(m, Gradient = ., expected = TRUE)# weights: 3 (2 variable)
initial value 162.196440
final value 133.205769
converged
predY <- pred %>%
do.call(rbind,.) %>%
as.data.frame() %>%
mutate(dominance=rep(gradient,1000)) %>%
pivot_longer(-c(dominance), names_to = "genome", values_to = "value")10.1.1 Plot responses to dominance
The estimated responses of genomes exhibiting significant positive and negative association with dominance.
#Get phylum colors from the EHI standard
phylum_colors <- genome_metadata %>%
left_join(read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv"), by=join_by(phylum == phylum)) %>%
arrange(match(genome, genome_tree$tip.label)) %>%
select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
slice(2:5) %>%
select(colors) %>%
pull()
getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="dominance", trend!="Neutral") %>%
select(genome,trend) %>%
left_join(predY, by=join_by(genome==genome)) %>%
filter(trend != "Neutral") %>%
#filter(genome %in% predY_asc) %>% #only display genomes with contrasting dynamics across treatments
group_by(genome, trend, dominance) %>%
summarize(value = mean(value, na.rm = TRUE)) %>%
left_join(genome_metadata, by=join_by(genome == genome)) %>%
ggplot(aes(x=dominance, y=value, group=genome, color=phylum)) +
geom_line() +
scale_color_manual(values=phylum_colors) +
facet_grid(fct_rev(trend) ~ phylum) +
labs(y="Genome abundance (log)",x="Dominance") +
theme(legend.position = "none") +
theme_minimal() +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 0.8,),
axis.line.x = element_line(size = 0.5, linetype = "solid", colour = "black"),
)10.2 Plot traits linked to trends
10.2.1 Test functional trait differences
Using Wilcoxon test with Benjamini-Hochberg procedure for FDR, identify functional traits that differ between genomes that are either positively or negatively association with dominance.
postestimate_test <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="dominance", trend!="Neutral") %>%
select(genome,trend) %>%
left_join(function_table %>% rownames_to_column(var="genome"), by=join_by(genome==genome)) %>%
pivot_longer(-c(genome,trend), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ trend)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)There are no traits that differ significantly between positively and negatively associated genomes.
10.2.2 Plot functional trait differences (not useful, as there are no differences)
postestimate_support %>%
rownames_to_column(var="genome") %>%
filter(trend != "Neutral") %>%
left_join(function_table %>% rownames_to_column(var="genome"), by=join_by(genome==genome)) %>%
pivot_longer(-c(genome,trend), names_to = "trait", values_to = "value") %>%
filter(trait %in% postestimate_test$trait) %>% #filter significantly different functions
filter(trait != "S0201") %>% #remove structural trait
ggplot(.,aes(x=value, y=trend, fill=trend, color=trend)) +
#geom_boxplot(color="black") +
geom_violin() +
scale_color_manual(values=c("#cf8888","#a2cc87"))+
scale_fill_manual(values=c("#cf8888","#a2cc87"))+
facet_grid(trait ~ ., space="free")+
theme_classic() +
labs(x="Metabolic Capacity Index",y="Metabolic function", fill="Trend", color="Trend") +
theme(strip.text.y = element_text(angle = 0),
strip.background = element_blank(),
axis.text.y=element_blank())10.2.3 Community functional predictions
function_table2 <- function_table %>% to.functions(., GIFT_db)
community_gifts <- predY %>%
group_by(dominance, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
select(-row_id) %>%
column_to_rownames(var = "dominance") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(function_table2,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="dominance")
})
community_gifts %>%
bind_rows() %>%
pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
mutate(dominance=as.numeric(dominance)) %>%
filter(trait=="B0218") %>%
ggplot(.,aes(x=dominance, y=value)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ splines::bs(x, 3), se = TRUE) +
theme_classic() +
labs(x="Dominance",y="Metabolic Capacity Index")
community_gifts %>%
bind_rows() %>%
pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
mutate(dominance=as.numeric(dominance)) %>%
#filter(trait=="B03") %>%
ggplot(.,aes(x=dominance, y=value)) +
geom_point() +
geom_quantile(quantiles = 0.05, formula = y ~ splines::bs(x, 3)) +
geom_quantile(quantiles = 0.5, formula = y ~ splines::bs(x, 3)) +
geom_quantile(quantiles = 0.95, formula = y ~ splines::bs(x, 3)) +
facet_grid(. ~ trait, ) +
theme_classic() +
labs(x="Dominance",y="Metabolic Capacity Index")